#####
# Replication for: "Can Political Speech Foster Tolerance of Immigrants?" by Schleiter, Tavits, and Ward.
# Figure S.6
#####

library(here)
library(data.table)
library(ggplot2)

# source the custom functions
source("functions.R")

# load the pooled data if not already in workspace
if(!exists("pooled")){
  pooled <- fread(here("data", "pooled.csv"))
}

# create the trichotomous education variable
pooled[ , edu3 := factor(
  fcase(
    edu2 %in% c("Less than high school", "High school graduate/ GED"), "Low",
    edu2 %in% c("2 year degree", "Some college"), "Medium",
    edu2 %in% c("4 year degree", "Masters degree", "Professional Degree (JD, MD)", "Doctorate"), "High"
  ), levels = c("Low", "Medium", "High")
)]


# Helper functions --------------------------------------------------------

# These three functions do the starch of the analysis. The first fits all 9 models, the second extracts the 9 marginal effects per model. The third does the battery of p-values
allCovEdufits <- 
  function(y, samp, data = pooled, chatty = T){
    if(samp != "Replication"){
      form <- as.formula(paste0(y, "~ CH + NM + CS + CH*edu3 + NM*edu3 + CS*edu3 + female + age + I(age^2) + race2 + pid2 + news2 + region2 + factor(wave)"))
    } else {
      form <- as.formula(paste0(y, "~ CH + NM + CS + CH*edu3 + NM*edu3 + CS*edu3 + female + age + I(age^2) + race2 + pid2 + news2 + region2"))
    }
    
    sample_use <- c("Main" = "replication == 0", "Replication" = "replication == 1", "Pooled" = "")[samp]
    
    if(samp != "Pooled"){
      mod <- lm(formula = form, data = data[eval(parse(text = sample_use))])
    } else {
      mod <- lm(formula = form, data = data)
    }
    
    if(chatty) print(summary(mod))
    return(mod)
  }

edu_ME_extract <- 
  function(mod, level = .95){
    out <- expand.grid("treatment" = c("CH", "NM", "CS"), "edu" = c("L", "M", "H"), stringsAsFactors = F); setDT(out)
    
    out[edu == "L" , c("est", "lower", "upper") := .(
      coef(mod)[c("CH", "NM", "CS")],
      confint(mod, level = level)[c("CH", "NM", "CS"), 1],
      confint(mod, level = level)[c("CH", "NM", "CS"), 2]
    )]
    
    for(n in 4:9){
      vars <- c(out$treatment[n], paste0(out$treatment[n], ":edu3", ifelse(out$edu[n] == "M", "Medium", "High")))
      est <- sum(coef(mod)[vars])
      se <- sqrt(vcov(mod)[vars[1], vars[1]] + vcov(mod)[vars[2], vars[2]] + 2*vcov(mod)[vars[1], vars[2]])
      set(out, i = n, j = "est", value = est)
      set(out, i = n, j = "lower", value = est - 1.96*se)
      set(out, i = n, j = "upper", value = est + 1.96*se)
    }
    
    return(out)
  }


# Fitting -----------------------------------------------------------------

edu_fits <- expand.grid(
  "y" = c("immPCA", "imm_neighbors2", "imm_increase2"),
  "samp" = c("Main", "Replication", "Pooled")
)
setDT(edu_fits)

edu_fits[ , fit := Map(allCovEdufits, y, samp, chatty = F)]
edu_me<- edu_fits[ , edu_ME_extract(mod = fit[[1]]), by = .(y, samp)]


# Plot prep ---------------------------------------------------------------


edu_me[ , yval := c("CH" = 1L, "CS" = 2L, "NM" = 3L)[treatment]]
edu_me[ , edu2 := factor(
  fcase(
    edu == "L", "Low",
    edu == "H", "High",
    edu == "M", "Medium"
  ),
  levels = c("Low", "Medium", "High")
)]
edu_me[ , y2 := factor(y, labels = c("Imm. Index (SD = 1)", "Imm. Neighbors (1-7)", "Increase Imm. (1-7)"))]
edu_me[ , samp2 := factor(samp, labels = c("Main Sample", "Replication Sample", "Pooled Sample"))]


# Plotting ----------------------------------------------------------------

(edu_plot <- ggplot(edu_me, aes(x = est, xmax = upper, xmin = lower, y = yval, color = edu2, shape = edu2)) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggstance::geom_pointrangeh(position = ggstance::position_dodge2v(.7, reverse = T)) +
    facet_grid(y2 ~ samp2) +
    scale_y_continuous(breaks = c(1,2,3), labels = c("Common\nHumanity", "Countering\nStereotypes", "Norms"), limits = c(.5, 3.5), expand = expansion(0,0)) +
    ylab(NULL) +
    xlab("Estimated Treatment Effect") +
    scale_color_grey(name = "Education", end = 0.7) +
    scale_shape(name = "Education") +
    theme_bw() +
    theme(legend.position = "bottom",
          axis.title.x = element_text(size = 10),
          axis.text = element_text(size = 10),
          strip.text = element_text(size = 10),
          legend.text = element_text(size = 10),
          legend.title = element_text(size = 10))
)
SI_width = 6.50127  

ggsave(here("SI", "figure_s6.pdf"), edu_plot, width = SI_width, height = SI_width*1.25)
